'''
结果：2.65236725 2.61366605 1.16279663 0.52913581 5.43101431 0.53439479 0.23676425 0.63855159 0.40290452
'''
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import least_squares
import pandas as pd


# 1. 读取Excel文件
def load_excel_data(file_path='attachment2_data.csv', sheet_name='Sheet1'):
    """读取Excel数据并返回DataFrame"""
    df = pd.read_excel(file_path, sheet_name=sheet_name)
    print(f"成功读取数据，维度：{df.shape}")
    return df


# 2. 数据清洗与格式转换
def preprocess_data(df):
    """数据预处理"""
    # 列名标准化（中英文兼容）
    column_map = {
        'temperature(℃)': 'temperature',
        'humidity': 'humidity',
        'solid content': 'solid_content',
        'proportion of hole area': 'hole_area'
    }
    df.rename(columns=column_map, inplace=True)

    # 处理缺失值（删除含NaN的行）
    df.dropna(subset=['temperature', 'humidity', 'solid_content', 'hole_area'], inplace=True)

    # 转换数据类型
    df = df.astype({
        'temperature': 'float',
        'humidity': 'float',
        'solid_content': 'float',
        'hole_area': 'float'
    })

    # 转换为NumPy数组（保留原始顺序）
    data_array = df[['temperature', 'humidity', 'solid_content', 'hole_area']].to_numpy()

    print(f"预处理后数据维度：{data_array.shape}")
    return data_array


# 读取数据（示例路径）
raw_data = load_excel_data(file_path='attachment2_data.csv')

# 预处理
processed_data = preprocess_data(raw_data)

# 在模型中调用
data = processed_data


def model_prediction(params, T, H, SC):
    """
    预测给定制备条件下的孔面积占比

    参数:
    params : list [k1, kH, k3, k4, k5, k6, S_C, S_S, alpha]
    T : float - 温度(℃)
    H : float - 湿度(%)
    SC : float - 固含量(%)

    返回:
    P_pred : float - 预测孔面积占比(%)
    """
    # =============================
    # 1. 参数解包与单位转换
    # =============================
    k1, kH, k3, k4, k5, k6, S_C, S_S, alpha = params
    T_k = T + 273.15  # 摄氏温度转开尔文

    # =============================
    # 2. 常量定义 (基于文献)
    # =============================
    # 初始质量 (kg)
    m_D0 = 24e-3
    m_S = 6e-3
    m_C = (m_D0 + m_S) * (SC / 100)  # 固含量转换为小数

    # 密度 (kg/m³)
    ro_D = 948
    ro_S = 1261
    ro_C = 1300

    # 物性参数
    A = 6.09451
    B = 2725.96
    C = 28.209
    r = 0.3413e-9  # 分子半径 (m)
    k_boltz = 1.380649e-23  # 玻尔兹曼常数

    # =============================
    # 3. 蒸发终止条件计算
    # =============================
    # 计算蒸发速率系数A0
    P_sat = 10 ** (A - B / (T_k + C)) * 0.133322  # Antoine方程计算饱和蒸气压
    humidity_factor = (1 - kH * H / 100)
    V = m_D0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积(m³)
    A0 = -k1 * P_sat * humidity_factor / (ro_D * V)

    # 求解蒸发停止时间 (DMF残留α%)
    def m_D(t):
        return m_D0 * np.exp(A0 * t)

    # 数值求解 t_stop: m_D(t) = α * m_D0
    t_stop = np.log(alpha) / A0

    # =============================
    # 4. 析出物计算
    # =============================
    def m_S_out(t):
        dissolved = min(m_S, m_D(t) * S_S)
        return m_S - dissolved

    def m_C_out(t):
        dissolved = min(m_C, m_D(t) * S_C)
        return m_C - dissolved

    def phi_C_out(t):
        return m_C_out(t) / (ro_C * V)

    # =============================
    # 5. 液滴动力学模型
    # =============================
    # DMF粘度计算 (不同温度)
    viscosity_table = {30: 0.8007, 40: 0.6560, 50: 0.5494}
    eta0 = viscosity_table.get(round(T), 0.6560) * 1e-6  # mPa·s → kPa·s

    def eta(t):
        return eta0 * (1 + 2.5 * phi_C_out(t))

    def D(t):
        return k_boltz * T_k / (6 * np.pi * r * eta(t))

    def v(t):
        return k3 * np.sqrt(D(t))

    # =============================
    # 6. 液滴碰撞微分方程
    # =============================
    def n_density(t):
        return m_S_out(t) / (ro_S * 4 / 3 * np.pi * r ** 3) / V

    def Z(t):
        return np.sqrt(2) * n_density(t) * np.pi * r ** 2 * v(t)

    def droplet_ode(m_s, t):
        """大液滴质量变化的微分方程"""
        dmdt = k4 * Z(t) + k5 * v(t) * (m_S_out(t) - m_s) / V * m_s
        return dmdt

    # =============================
    # 7. 求解微分方程
    # =============================
    t_eval = np.linspace(0, t_stop, 50)  # 50个时间点
    sol = odeint(droplet_ode, 0, t_eval)  # 初始质量=0
    m_s_final = sol[-1][0]  # 最终大液滴质量

    # =============================
    # 8. 孔面积占比预测
    # =============================
    droplet_volume = m_s_final / ro_S
    P_pred = k6 * (droplet_volume / V) * 100  # 转换为百分比
    print(f"孔面积占比预测：{P_pred}")
    return P_pred

def residual(params):
    """ 计算所有数据点的残差 """
    residuals = []
    for T, H, SC, P_exp in data:
        P_pred = model_prediction(params, T, H, SC)
        residuals.append(P_pred - P_exp)
    return np.array(residuals)

# 参数初值及边界
initial_guess = [1.0, 0.9, 1.0, 1.0, 1.0, 1.0, 0.26, 0.26, 0.3]  # 基于文献的合理初值
bounds = ([0.1, 0.001, 0.1, 0.1, 0.1, 0.1, 0.01, 0.01, 0.0],
          [10.0, 1.0, 10.0, 10.0, 10.0, 10.0, 1.0, 1.0, 0.41177578])

# 执行优化
result = least_squares(residual, initial_guess, bounds=bounds, method='trf')
optimized_params = result.x
print(optimized_params)
